† Corresponding author. E-mail:
Project supported by the National Natural Science Foundation of China (Grant No. 11547009), the Young Scientists Fund of the National Natural Science Foundation of China (Grant No. 11602062), the Natural Science Foundation of Guizhou Province, China (Grant No. 2012/2166), and the Research Foundation of Guizhou University for Talent Introduction, China (Grant No. 2011/02).
The anisotropy in the particle systems of different packing structures affects the sound velocity. The acoustic propagation process in four kinds of packing structures (denoted as S45, H60, S90, and D) of two-dimensional granular system is simulated by the discrete element method. The velocity vtof obtained by the time of flight method and the velocity vc obtained from the stiffness tensor of the system are compared. Different sound velocities reflect various packing structures and force distributions within the system. The compression wave velocities of H60 and S90 are nearly the same, and transmit faster than that of D packing structure, while the sound velocity of S45 is the smallest. The shear wave velocities of S45 and H60 are nearly the same, and transmit faster than that of D packing structure. The compression wave velocity is sensitive to the volume fraction of the structure, however, the shear wave velocity is more sensitive to the geometrical structure itself. As the normal stress p is larger than 1 MPa, vtof and vc are almost equal, and the stiffness tensors of various structures explain the difference of sound velocities. When the normal stress is less than 1 MPa, with the coordination number unchanged, the law vtof ∝ p1/4 still exists. This demonstrates that apart from different power laws between force and deformation as well as the change of the coordination number under different stresses, there are other complicated causes of vtof ∝ p1/4, and an explanation of the deviation from vtof ∝ p1/6 is given from the perspective of dissipation.
In granular assemblies, there are inhomogeneous contact forces and geometrical structures on the scale of granular particles, which makes the granular system show complicated properties such as nonlinearity, disorder, and anisotropy,[1] and the macroscopic description of the material behavior is very difficult. The inhomogeneous stress distribution in quasistatic granular packings will cause the contact force between some particles to be greater than the average stress; these particles usually exist in the form of chains, and support the whole system.[2] It is still not clear that how the heterogeneity of the stress distribution and the existence of the force chain affect the properties of the granular system, whether or not the granular system can be regarded as an elastomer. If it can be regarded as an elastomer, it is also controversial whether the elastic theory is applicable to describe the granular system at various scales.[2] To obtain the microstructure evolution and macroscopic mechanical properties of a system, the most intuitive research methods are the photoelastic experiment,[3] x-ray,[4] and nuclear magnetic imaging,[5] which are expensive. Moreover, these methods also have some difficulties to detect the real granular system. Besides these methods, acoustic propagation experiments can be used to understand the relationships among the micro-stress distributions, geometrical structures, and macroscopic properties. The nonlinearity and disorder of the granular system play important roles in the propagation of mechanical energy, and the microstructure and stress distribution correspond to different propagation modes.[6] The acoustic propagation experiment provides an effective way to detect the internal structure of the real granular system. For example, the low frequency acoustic detection experiment can determine the nonlinear elastic property and structural anisotropy of the granular media,[7–9] while short-wave multiply scattered sound waves enable one to detect tiny changes of the contact network configuration at the grain level,[10,11] as well as the use of acoustic detection techniques to study seismic triggering[12] and the formation of shear bands, etc.[13]
In the acoustic experiment, the determination of the sound velocity is usually done by the time of flight method.[14,15] This method obtains the time and the corresponding velocity of the sound wave in the granular system through the real path because sound waves favor speed and amplitude along force chains.[16,17] The waveform and velocity obtained by the time of flight method reflect information such as stress, strain, stiffness tensor, and attenuation along the wave propagation path. There is another method to obtain the sound velocity apart from the time of flight method. If the granular system is regarded as a continuous medium, the elastic potential energy of the system can be calculated according to the contact relation based on the effective medium theory (EMT)[8] or granular solid hydrodynamics (GSH) theory.[15] Then, the stress can be obtained by taking partial derivative of the potential energy with respect to strain, and the sound velocity can be obtained from the fourth-order stiffness tensor.[18,19] EMT and GSH theories usually deal with homogeneous cases without considering the inhomogeneity of geometry and stress distribution, and the given stiffness tensor and sound velocity are the average effect of the whole system.[15]
Whether from the time of flight method or the stiffness tensor, the sound velocity is the response of the system stress structure. For example, the relationship between the normal stress p and sound velocity of the disordered packing granular system obtained by the time of flight method vtof meets vtof ∝ p1/6 while this dependence is at higher pressures, but the behavior at lower confining pressures is more like vtof ∝ p1/4.[7,15,20] However, the sound velocity calculated from the stiffness tensor vc meets the relation vc ∝ p1/6.[2,8,21,22] The compression waves are more sensitive to the stress-induced anisotropy, whereas the shear waves are more sensitive to the fabric anisotropy.[9] The two-dimensional photoelastic experiment shows that the sound velocity v has a nonlinear behavior with the applied force F as v = KFα, where the exponent α depends on the contact number of the grains, while the prefactor K brings information about the anisotropy of the medium.[17,23] Therefore, the study of the relationship between sound velocity and the internal structure of particles can not only detect a change in the structure of the granular system, but also recognize its nature, and it is of great interest when trying to predict abrupt events, namely, avalanches, earthquake, or unjamming event.[16,24]
The anisotropic granular system[1] is obtained experimentally with different preparation methods[9] (preparation protocols of rain deposition, decompaction) or different stress loading methods; however, it is difficult to measure the internal stress distribution of the system. Usually, the force on the boundary wall is used to substitute the stress distribution in the system,[9] and the calculation of stiffness tensor has also been carried out by some simplification and approximation.[14] The limitation of experimental means makes it difficult to accurately describe the particle system. The acoustic experiment, especially the acquisition of the stable shear wave, is also more difficult and higher requirement for equipment is needed.
The discrete element method (DEM) can obtain the force and deformation of each particle, and simulate the sound wave through the simple harmonic vibration or simple harmonic force of the boundary wall, making it convenient to obtain the waveform and sound velocity. Furthermore, DEM can obtain the required physical quantities of EMT and GSH theories to explain acoustical properties of granular systems, for example, the coordination number used to explain the relation v ∝ p1/4,[25] the volume and shear modulus used to calculate the sound velocity,[8] and so on.
In this paper, four different packing structures (S45, H60, S90, and D) in the two-dimensional photoelastic experiment are simulated by DEM.[17] The relationship between compressional and shear wave velocities and normal stress is calculated by the time of flight method. The influence of different packing structures on the sound velocity under the same normal stress is analyzed. The stress tensor and the stiffness tensor under different packing structures are calculated. The wave velocities obtained by the stiffness tensor are compared with those by the time of flight method. The variation of sound velocity with stress and packing structures can be explained by the difference of system stiffness and stress tensor. The results are consistent with the existing experiments and theories.[17] The change of sound velocity reflects the difference of stress distribution and geometrical structure in the system, which can provide a reference for acoustically probing the structure of system.
Four kinds of packing structures used to calculate the sound velocity are shown as S45, H60, S90, and D (disordered) in Fig.
The contact between particles is simulated by the simplified Hertz–Mindlin–Deresiewicz model. The force and torque expressions are listed in Table
![]() | Table 1. Force and torque between particles. . |
The parameters of the simulations are listed in Table
![]() | Table 2. Parameters of the simulations. . |
Using walls as boundary conditions, the top wall is subjected to normal stress, while the other three walls are fixed in the loading process. When the system achieves equilibrium completely under a certain normal stress, the sound wave is simulated by a simple harmonic vibration of the bottom wall. The simple harmonic vibrations in the vertical and horizontal directions represent the compression and shear waves, respectively, as shown in Fig.
![]() | Fig. 3. (color online) Input and output waveforms: (a) input waveform, (b) output compression waveform, and (c) output shear waveform. |
Taking the time at 1/10 of the first wave peak as the flight time of the sound wave, the distance between the top and bottom walls is divided by the time to give the velocity. As shown in Fig.
The compression and shear wave velocities of four packing structures under the normal stresses from 200 kPa to 3 MPa are shown in Figs.
![]() | Fig. 4. (color online) Compression wave velocities under different normal stresses and for different structures: (a) S45, (b) H60, (c) S90, (d) D. |
There are many theories and methods to convert the microscopic contact properties between particles into the macro-mechanical properties, such as stress, strain, stiffness tensor, and so on, which are used to deal with the wave propagation in the classical theory.[1,26–28] As with the sound velocity obtained by the time of flight method, the average stress and stiffness tensor from the microscopic statistics reflect the internal stress distribution and geometric structure of the system. The influence of stress and structural inhomogeneity on the acoustic propagation can be obtained by calculating the stress tensor and the stiffness tensor.
At the contact point, the physical quantity
![](cpb_27_10_104501/cpb_27_10_104501_ieqn10.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn11.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn12.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn13.gif)
The relationship between velocity and stiffness tensor[14,29] is
The velocities by the time of flight method and stiffness tensor are related to the packing structures. The compression wave velocities of H60 and S90 are larger than those of S45 and disordered structure under the same normal stress, and the relative velocity difference can reach 30%. For instance, at the normal stress of 3 MPa,
The relationship between shear wave velocities and packing structures are different. The shear wave velocity of S45 is nearly the same as that of H60, and bigger than that of disordered packing. For example, at the normal stress of 50 kPa,
The velocity difference between various structures can also be explained by different stiffness tensors. Different structures mean different geometries and force distributions, and then determine the stiffness tensors. From the ideal of a continuum, velocity depends on stiffness, and in fact, the velocity difference is the stiffness difference determined by the geometry. The change in stiffness determines the change in velocity, especially for greater normal stress (more than 1 MPa).
In addition to stiffness, the difference of velocity can be interpreted by the dispersion of elastic waves in granular phononic crystals. In the case of hexagonal close packed structure, as shown in Fig.
![](cpb_27_10_104501/cpb_27_10_104501_ieqn22.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn23.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn24.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn25.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn26.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn27.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn28.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn29.gif)
![](cpb_27_10_104501/cpb_27_10_104501_ieqn30.gif)
![]() | Fig. 6. (color online) Contact model of hexagonal close packed structure: (a) normal spring, (b) tangential spring, and (c) rolling spring. |
The dispersion relation gives the velocity of the shear wave and that of the compression wave, in the long-wave limit:
![](cpb_27_10_104501/cpb_27_10_104501_ieqn31.gif)
For the square phononic crystal (S90), the compression wave velocity is[30]
The ratio of compression wave velocities for H60 and S90 is
For S45, because of the symmetry, the first two equations of motions are the same, and the velocities of the compression and shear waves should be equal, which is in agreement with Figs.
The velocities from the lattice dynamics are ideally given, and the system is always symmetric and homogeneous. However, for a realistic granular system, it is difficult to prepare a sample in which the force distribution is the same as its geometry even by simulation. That is the reason for the deviation between the velocities by the time of flight method and those by the stiffness or lattice dynamics.
The relationship between the normal stress and velocity derived from the stiffness tensor is vc ∝ p1/6 as shown in Figs.
The shear wave velocities by the time of flight method have an overall deviation from the velocities (Fig.
![]() | Fig. 7. (color online) Force distributions in samples of different preparation methods: (a) Fx in sample 1, (b) Fx in sample 2, (c) Fy in sample 1, (d) Fy in sample 2. |
Figures
![]() | Fig. 8. (color online) Influence of damping on the shear velocity: (a) Cn = 0.75, Ct = 0.5, and (b) Cn = 0.5, Ct = 0.25. |
The 1/6 exponent of Hertz scaling is limited to homogeneous and non-dissipative samples, however, dissipation is the nature of reality in the particle system. When the particle system achieves balance under external force, the strains and the relative positions between particles are unchanged. In the state of equilibrium, the system can be considered as an example of elastomer. When the system is disturbed by waves, the relative velocities and positions between the particles emerge. As a result of the relative velocity and position, the force chains in the system have the possibility of relaxation, then the dissipation emerges, which causes the system to deviate from elastomer, that is, deviate from the 1/6 exponent of Hertz scaling.
The fluctuation of energy under different normal stresses during the shear wave propagation is shown in Fig.
![]() | Fig. 9. (color online) Fluctuation of energy during wave propagation: (a) kinetic energy of particles, and (b) ratio of kinetic and potential energy. |
For small normal stress p = 100 kPa, the kinetic energy of particles is larger than the input wave energy, which means the difference of the energy comes from the initial elastic potential energy of the system. That is to say, some of the elastic potential energy stored in force chains is dissipated by the wave disturbance, especially in the tangential direction, and the number of the broken tangential springs at every time step is shown in Fig.
In Fig.
From another point of view, in the contact model used in this paper, cn and ct are used to balance the relative velocity and position between the particles. The larger the values of cn and ct, the more difficultly the particle moves; in other words, for the static system, cn and ct contribute to prevent the dissipation caused by the relative motions of the particles, that is why the shear wave velocity tends to keep accordance with the 1/6 exponent of Hertz scaling as cn and ct mounted, as shown in Figs.
In the end, velocities of four packing structures reflect different geometries and force distributions; the compression wave velocity is sensitive to the volume fraction of the structure, and the shear wave velocity is more sensitive to the geometrical structure itself. For the stress dependence of the sound wave velocity, the dissipation is the probable cause of the deviation from 1/6 exponent of Hertz scaling, and further research will be done in the future.
[1] | |
[2] | |
[3] | |
[4] | |
[5] | |
[6] | |
[7] | |
[8] | |
[9] | |
[10] | |
[11] | |
[12] | |
[13] | |
[14] | |
[15] | |
[16] | |
[17] | |
[18] | |
[19] | |
[20] | |
[21] | |
[22] | |
[23] | |
[24] | |
[25] | |
[26] | |
[27] | |
[28] | |
[29] | |
[30] | |
[31] | |
[32] | |
[33] |